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Abstract. The complex wave representation (CWR) converts unsigned 2D distance transforms 
into their corresponding wave functions. Here, the distance transform S(X) appears as the phase of 
the wave function 4>{X) — specifically, c/>(X) = exp ( i§12Q. j where r is a free parameter. In this work, 
we prove a novel result using the higher-order stationary phase approximation: we show convergence 
of the normalized power spectrum (squared magnitude of the Fourier transform) of the wave function 
to the density function of the distance transform gradients as the free parameter r — ¥ 0. In colloquial 
terms, spatial frequencies are gradient histogram bins. Since distance transform gradients carry only 
orientation information (as their magnitudes are identically equal to one almost everywhere), the 2D 
Fourier transform values mainly lie on the unit circle in the spatial frequency domain as r — > 0. The 
proof of the result involves standard integration techniques and requires proper ordering of limits. 
Our mathematical relation indicates that the CWR of distance transforms is an intriguing, new 
representation. 
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1. Introduction. Euclidean distance functions (more popularly referred to as 
distance transforms) are widely used in many domains |17U22| . An important subset — 
point-set based distance functions — also finds application in many domains with com- 
puter vision being a prominent example |22[ 111! [20] . Since distance transforms allow 
us to transition from shapes to a scalar field, problems such as shape registration are 
often couched in terms of rigid, affine or nonrigid alignment of distance transform 
fields, where the shapes are parameterized as a set of points [IS]. In medical imaging, 
they are used in the construction of neuroanatomical shape complex atlases based on 
an information geometry framework [2]. 

Even when one begins with a set of closed curves (as a shape template for ex- 
ample), the curves are often discretized to yield a point-set prior to the application 
of fast sweeping |25| and other distance transform estimation methods. Signed and 
unsigned distance transforms are deployed in 3D as well with their zero level-sets cor- 
responding to surfaces. Furthermore, medial axis methods and skeletonization often 
involve distance transform representations |13| . In the domain of computer vision, the 
gradient density function is popularly known as the histogram of oriented gradients 
(HOG). Since the advent of HOG a few years ago, gradient density estimation has 
risen in prominence and is employed in human recognition systems [¥]. 

The distance transform for a set of K discrete points Y = {Yk E R D },k £ 
{1, . . . , K} where D is the dimensionality of the point-set is defined as 

(1.1) S(X) = mm\\X-Y k \\, 

k 

where X G ft is a closed bounded domain in R 15 . In this article, we are only concerned 
with D = 2. 
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In computational geometry, Euclidean distance functions correspond to the Voronoi 
problem [5] and the solution S(X) can be visualized as a set of cones (with the cen- 
ters being the point-set locations {Yfe}). The distance transform satisfies the static, 
non-linear Hamilton- Jacobi equation 



almost everywhere, barring the point-set locations and the Voronoi boundaries where 
it is not differentiable [T7 | fT8 l l2~T]. Here VS = (S x ,S y ) denotes the gradients of S 
and || • || represents its Euclidean magnitude. Furthermore S(X) = at the point-set 
locations. Following the wave optics literature, one can envisage light waves simulta- 
neously emanating from the given point sources and propagating with a velocity of 
one in all directions. The value of S at a grid point Xo, namely S(Xq), corresponds 
to the time taken by the first light wave (out of the K light waves) to reach the grid 
location X . Driven by this optics analogy, when we express 5* as the phase of a wave 
function (f> as in 



we made an intriguing empirical observation. The power spectrum of the wave func- 
tion approximates the density function of the gradients of the distance transform as 
the parameter r in Equation 11.31 tends to zero. In this paper, we formally prove this 
result. We refer to this wave function <j> which satisfies the phase relation with S as 
the Complex Wave Representation (CWR) of distance transforms. 

2. Main Contribution. The centerpiece of this work is to provide a useful 
application of the stationary phase method, wherein we show an equivalence between 
the density function of the gradients of the distance function VS = (S x , S y ) and the 
power spectrum (squared magnitude of the Fourier transform) of the CWR ((f)) as 
the free parameter t (in Equation II. 311 approaches zero. Here, the density function 
of the gradients is obtained via a random variable transformation of a uniformly 
distributed random variable W (over the bounded domain £1) using the gradients 
VS = (S x , S y ) as the transformation functions. In other words, if we define a random 
variable Z = WS(W) where the random variable W has a uniform distribution on 
a closed bounded domain £1 C M 2 , the density function of Z represents the density 
function of the gradients of the distance transform. 

As the norm of the gradients VS is defined to be 1 almost everywhere (from 
Equation ll.2p , we observe that the density function of the gradients is one-dimensional 
and defined over the space of orientations. Section[3]provides a closed-form expression 
for this density function. As the gradients are unit vectors, we notice that the Fourier 
transform values of the CWR ((f)) lie mainly on the unit circle and this behavior 
tightens as r — > 0. Specifically, if F T (r,u}) represents the Fourier transform of <f> in 
the polar coordinate system at a given value of t, Theorem 14.21 demonstrates that if 
f =/= 1, then lim T ->o F T (f, uS) = 0. 

Our main result is established in Theorem 15.11 where we show that the power 
spectrum of the wave function </> when polled close to the unit circle, is approximately 
equal to the density function of the distance transform gradients, with the approxi- 
mation becoming increasingly exact as r — > 0. In other words, if P(u) denotes the 
closed-form density of the gradients defined over the orientation uj and if P T (f,uS) 
corresponds to the power spectrum of (f> represented in the polar coordinate system 



(1.2) 



||VS|| = 1 



(1.3) 




DISTANCE TRANSFORM DENSITY ESTIMATION 



3 



at a given value of t, Theorem 15. II constitutes the following relation 



(2.1) limlim / < P r (r,uS)rdr \ duj = / P(oj)duj 

' U)Q I J 1 — 5 I J UJQ 



<5->0 t->0 



for any (small) value of the interval measure A on w. We show this result using the 
higher-order stationary phase approximation, a well known technique in asymptotic 
analysis [23] . Through the pioneering works of Jones and Kline [T3j , Olver [TS] , Wong 
|23j . McClure and Wong Q3], among others, the stationary phase approximation has 
become a widely deployed tool in the approximation of oscillatory integrals. Our 
work showcases a novel application of the stationary phase method for estimating 
the probability density function of distance transform gradients. The significance of 
our mathematical result is that spatial frequencies become histogram bins and hence 
the power spectrum P T can serve as a gradient density estimator at small, non-zero 
values of r. We would like to emphasize that our work is fundamentally different 
from estimating the gradients of a density function [8 and should not be semantically 
confused with it. 

2.1. Motivation from quantum mechanics. Our new mathematical rela- 
tionship is motivated by the classical-quantum relation, wherein classical physics is 
expressed as a limiting case of quantum mechanics [TU1 |B]. When S is treated as the 
Hamilton- Jacobi scalar field, the gradients of S correspond to the classical momentum 
of a particle [9]. In the parlance of quantum mechanics, the squared magnitude of 
the wave function expressed either in its position or momentum basis corresponds to 
its position or momentum density respectively. Since these representations (either 
in the position or momentum basis) are simply (suitably scaled) Fourier transforms 
of each other, the squared magnitude of the Fourier transform of the wave function 
expressed in its position basis is its quantum momentum density. However, the time 
independent Schrodinger wave function 4>{x, y) (expressed in its position basis) can 

be approximated by exp 

^iSl^vl^ as t — } [6J. Here r (treated MS M free parameter in 

our work) represents Planck's constant. Hence the squared magnitude of the Fourier 

transform exp ( ^ s ^ v ^> ^ corresponds to the quantum momentum density of S. The 

principal results proved in the article (Theorem 15.11 and Proposition 16. 5[) state that 
the classical momentum density (denoted by P) can be expressed as a limiting case 
(as t — y 0) of its corresponding quantum momentum density (denoted by P T ), in 
agreement with the correspondence principle. 

3. The Distance Transform Gradient Density Function. As mentioned 
above, the geometry of the distance transform corresponds to a set of intersecting 
cones with the origins at the Voronoi centers [5] . The gradients of the distance trans- 
form (which exist globally except at the cone intersections and origins) are unit vectors 
and satisfy Equation 11.21 Therefore the gradient density function is one-dimensional 
and defined over the space of orientations. The orientations are constant and unique 
along each ray of each cone. Its probability distribution function is given by 

(3.1) F{6 < 6 < + A) = - f [ dxdy 



<e+A 



where L is the area of the bounded domain Q. We have expressed the orientation 
random variable as O = arctan (^^j ■ The probability distribution function also 
induces a closed-form expression for its density function as shown below. 
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Let fl C M 2 denote a polygonal grid such that its boundary dfl is composed of a 
finite sequence of straight line segments. The reason for restricting only to polygonal 
domains with boundaries made of line segments will become clear when we discuss 
Theorem [OJ Let the set Y = {Y k G R 2 ,fc G {1, ...,K}} be the given point-set 
locations and let Yk — (xk,Vk)- Then the Euclidean distance transform at a point 
X = (x, y) G fl is given by 

(3.2) S(X) = min \\X-Y k \\= mm(^/(x - x k ) 2 + (y - y k ) 2 ). 

k k 

Let T>k, centered at Yk, denote the k th Voronoi region corresponding to the input 
point Yk. T>k can be represented by a Cartesian product [0, 2w) x [0, Rk(&)] where 
Rk(0) is the length of the ray of the fc th cone at an orientation 6. If a grid point 
X = (x, y) G (Yfc + T>k), then S(X) = \\X — Each 2\. is a convex polygon whose 
boundary dT>k is also composed of a finite sequence of straight line segments as shown 
in Figure [XU 



Figure 3.1. Voronoi diagram of the given K points. Each Voronoi boundary is composed of 
straight line segments. 



Note that even for points that lie on the Voronoi boundary where the radial length 
equals R^ (9) , the distance transform is well defined. The area L of the polygonal grid 
f2 is given by 



k=1 Jo Jo k=1 Jo 



(3.3) L =? A rdrdO = > / -^d0 



With the above set-up in place, after recognizing the cone geometry at each Voronoi 
center Yk, Equation 13. II can be simplified as 

1 K r 6+A r R k (8) i K r e + A r>2(n\ 

(3.4) F{6 < 6 < 6 + A) = - V / / rdrdO = - V / ^^dO. 

L k=1 Je Jo L k =i Je 2 

Following this drastic simplification, we can write the closed-form expression for the 
density function of the unit vector distance transform gradients as 

(3-5) P(0) = Inn ^ ^ + A ) = I f 

y J w a^o A L ^ 2 

fc=i 
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Based on the expression for L in Equation 13.31 it is easy to see that 

(3.6) / P(B)d0 = 1. 

Jo 

Since the Voronoi cells are convex polygons [5j, each cell contributes exactly one 
conical ray to the density function on orientation. 

4. Properties of the Fourier Transform of the CWR. Since the distance 
transform is not differcntiable at the point-set locations {Yk}^ =1 and also along the 
Voronoi boundaries dT>k, Vfe (a measure zero set in 2D), we restrict ourselves to the 
region which excludes both of them. To this end, let < e < \ be given. Let the region 

r>| centered at Yj. be represented by the Cartesian product [0, 27r) x [R^\9) 7 RP'(9)] 
where, 

R k 1) (9) = eR k {9) and 
(4.1) R k 2) (9) = (l-e)R k (9). 

The length of the ray at the orientation 6 in V k equals Rj£ (9) — R k (9). Note that 
in the definition of T>1 , we have explicitly removed the source point Yfe where the ray 
length r(9) = and the boundary of the Voronoi cell where r(9) = Rk(9) as shown 
in Figure |4~T1 




Figure 4.1. Region that excludes both the source point and the Voronoi boundary. 



Define the grid 

K 

(4.2) n*=\J(Y k +vi). 

fc=i 

Its area L e equals 

K ~2tt rR^He) K i-2-k P 2 

(4.3) m rdrd9 = (l-2e)J2 



&l d 9. 



From Equation 13.31 we get L e — (1 — 2e)L and hence lim^o L c = L. 
Let V = VJ7-. Define a function F :KxKxl->Cas 

(4.4) F'(u, „, r )s ^ll exp («fej4) exp 

CI' 

For a fixed value of t, define a function : R x K —> C as 

(4.5) F*(u,v) = F%u,v,t). 



(i 
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Note that F* is closely related to the Fourier transform of the CWR, = exp (— ) 
[Tj . The scale factor 2 * is the normalization factor such that the £ 2 norm of i 7 !^ is 
1 as seen in the following Lemma (with the proof given in Appendix [X| . 
Lemma 4.1. With defined as above, F* e L 2 (R 2 ) and \\F*\\ = 1. 

Consider the polar representation of the spatial frequencies (it, v) namely u = 
rcos(w) and v — rsin(o;) where f > 0. For (x, y) G (Y/. + ^ ^ x ~~ x k — rcos(9) 
and y — yk — rsin(6) where r g [11^(9), (6)]. Then Equation 14.41 can be rewritten 



K 



(4.6) F^(r,ij) = J2CkIk(r,i 

where 



k=l 



(4.7) C fc = exp | — ^- [f cos(o;)a;A; + f sin(w)y fe ] 
and 

i r 2v r R * )(e) (i i 

(4.8) I k (r,u)) = / / exn<-r[l-rcas(6-w)]}rdrd6. 

27TT/ e J J R w {e) [t J 

With the above set-up in place, we have the following theorem, namely, 
Theorem 4.2. [Circle Theorem] Iff^l, then, 

(4.9) Iim F*(f,u;) = 0, 



o 



for any < e < i 



2 ' 

4.1. An Intuitive Examination of Theorem 14.21 Before we furnish a rigor- 
ous proof for the aforementioned theorem, we provide an intuitive picture of why 
the statement is true. Observe that the first exponential exp fi^fcsl^ m Equa- 
tion |L4] is a varying complex "sinusoid" and the second exponential exp ^ i{ux+vy) \^ 
in Equation 14.41 is a fixed complex sinusoid at frequencies — and - along the x- 
and y-coordinate axes respectively. When we multiply these two complex exponen- 
tials, at low values of t, the two sinusoids are usually not "in sync" and cancella- 
tions occur in the integral. Exceptions to the cancellation happen at locations where 
VS = (S x ,S y ) — (u,v), as around these locations, the two sinusoids are in perfect 
sync. Since IIV^H = 1 for distance transforms, strong resonance occurs only when 
u 2 + v 2 = 1 (f = 1). When f ^ 1, the two sinusoids tend to cancel each other out as 
r — > 0, resulting in F^ becoming zero at these locations. 

4.2. Proof of Theorem 14.21 Having given an intuitive picture of why Theo- 
rem [L2] holds true, we now proceed with the formal proof. As each Ck is bounded, it 
suffices to show that if f ^ 1, then lim T _j.o Ik{f, to) = for all Ik- 
Proof. Consider the integral 

(4.10) I(r,w) = - — r / exp<^ -r [1 -fcos(0- w)] \rdrd6, 



^TTTl e J Q J R (1) {0 ) 1 t 
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where RW(6) = eR{6) &ndR^{9) = (l-e)R(6). Let the region [0, 2tt)x [R^(9), R<- 2 \9)} 
be denoted by T> t . R(9) is defined in such a way that the boundary of D e consists 
of a finite sequence of straight line segments as in the case of each T>t- Notice that 
V e doesn't contain the origin (0, 0). In order to prove Theorem 14.21 it is sufficient to 
show that lim T _j.o ^(^j w ) = 0. 

Let p(r, 9; r, uj) = r(l — f cos(6 — uj)) denote the phase term of I in Equation 14. 101 
for a given f and uj. The partial derivatives of p{r, 9] f, u>) (with f and uj held fixed) 
are given by 



(4.11) 



dp 
Or 



= 1 — f cos(9 — uj), 



dp 
86 



= rr sm( V — uj) 



Since T> e is bounded away from the origin (0, 0), Vp is well-defined and bounded and 
equals zero only when f = 1 and 9 = lj. Since r ^ 1 by assumption, no stationary 
point exists (Vp ^ 0) and hence we can expect I(r,ui) — > as r — > O [12l [24] . 
Below, we show this result more explicitly. 



Define a vector field u(r, 9; r, lj) 



_ Vp 



at a fixed value of f and uj. Note that 



V 

(4.12) 



u(r, 9; f, uj) exp 



ip(r, 9; f, uj) 



(V • u(r, 9; r, uj)) exp 



/ ip(r,9;r,oj) 



i (' ip(r,6:r,ujy 
— exp I | r 



where the gradient operator V 
we get 



a_ 

dr ' r 00 1 



Inserting Equation l4.12l in Equation l4.10l 



(4.13) 
where 

(4.14) I {2 \f,uj) 



I(f,uj) = I^\r,uj)-I^\P,uj), 



2nil e 



V • u(r, 9\ f, uj) exp 



ip{r,V;r,uj, 



drd9, and 



1 



, );/ , /; (V-n(r,9;r,uj))e Xp ( lpir ' ^ Uj) \ 



Consider the integral I^(f,ui). From the divergence theorem, we have 



(4.15) 



2nil e 



, , ip(r,6;f,u})\ , 
(u n) exp I I ds 



where T is the positively oriented boundary of 2? e , s is the arc length of T and n is the 
unit outward normal of T. The boundary T consists of two disjoint regions, one along 
r(9) = RW(6) and another along r(9) = RW{9). If the level curves of p(r, 9\ r, uj) 
are tangential to T only at a discrete set of locations giving rise to stationary points 
of the second kind [531 [Ml HI] — m other words, if p(r, 9; f, uj) is not constant along 
the boundary T for any contiguous interval of 9 — then, using the one dimensional 
stationary phase approximation [16] , I^'(f, uj) can be shown to be 0(^/t) and 
hence converges to zero as r — > 0. Since the boundary of T> e is composed of straight 
line segments (specifically not arc-like), we can show that the level curves of p(r, 9; r, uj) 
cannot overlap with T for a non-zero finite interval. (The next paragraph takes care 
of this technical issue and can be skipped without loss of continuity.) 
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The level curves of p(r, 9;f, ui) are given by R(9)(l — f cos(# — u)) — c, where c 
is a constant. Recall that each of the two disjoint regions of T is composed of a finite 
sequence of line segments. For the level curves of p(r, 9; f , uj) to coincide with T over 
a non-zero finite interval, y(9) = R(9) sin(0) = !_f cos(e-^) ano - x (&) ~ ^fi) cos ($) = 
i_f cos(e-ct)) should satisfy the line equation y = mx + b for some slope m and slope- 
intercept b, when 9 varies over some contiguous interval 9 £ [#i,#2]- Plugging in the 
value of y{9) and x{9) into the line equation and expanding cos(# — lj), we have 

(4.16) csin (9) = mccos(9) + b — bf[cos(9) cos(w) + sin(#) sin(cj)]. 
Combining the terms, we get 

(4.17) sin(0)[c + bf sin(o;)] — cos((9) [mc — bfcos(u>)] = b. 

By defining Ai = c + brsm(uj) and A2 = — (mc — 6f cos(w)), we see that sin(#) and 
cos(#) need to satisfy the linear relation 

(4.18) Ai sin(0) + A 2 cos(0) = b 

for 9 € [6*1,^2] in order for the level curves of p(r, 9; f, uS) to overlap with the piece- 
wise linear boundary T. As Equation 14.181 cannot be true for a finite interval of 9, 
JW(f, oS) = 0{y/r) as t — y and hence converges to zero in the limit. 

Now l( 2 \r,us) has a similar form as the original I(f,ui) in Equation 14.101 with 
r replaced by gt(r,9;r,w) = (V • u). Letting ui(r,d;f,u) = ^^ gi(r,9; f, cj), from 
Equation 14.121 and the divergence theorem, we get 



1 (r ' W) = 2^ y r (Ul n) CXP { r ) dS 



T 



ip(r, 



(4.19) + 2^FJJ ft 7 • Mr, 9;r,u)))exp ( r ) drdfl. 

v 

As I^ 2 \f, ui) = 0(t), it converges to zero as r — > 0. Applying the obtained results to 
Equation 2331 we see that I{f,bS) (and also lfe(f, lS) defined in Equation I4.8[) — > as 
t — V which completes the proof. 
□ 

Since Theorem 14. 2 1 is true for any < e < i, it also holds as e — > 0. As a corollary, 
we have the following result: 

Corollary 4.3. Iff^l, then 

(4.20) lim lim ,F T £ (f,w) = 0. 

e-s-0 r->-0 

5. Spatial Frequencies as Gradient Histogram Bins. We now show that 
the squared magnitude of the Fourier transform of the CWR (0) when polled close to 
the unit circle (f = 1) is approximately equal to the density function of the distance 
transform gradients (P) with the approximation becoming increasingly tight as r — > 0. 

The squared magnitude of the Fourier transform — also called its power spectrum 
PQ — is given by 



(5.1) 



P*(r,uj) = |F r e (f, W )| 2 = F T e (f, W )F r ^). 
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By definition, P*(r,uj) > 0. From Lemma l4.11 we have 

r>1n />0O 

(5.2) / / P*(r,u)rdrdu = 1 



independent of r. Hence, P^{f , w) can be treated as a density function for all values of 
r. We earlier observed that the gradient density function of the unit vector distance 
transform gradients is one-dimensional and defined over the space of orientations u. 
For P^(f, oj) to behave as an orientation density function, it needs to be integrated 
along the radial direction f . Since Theorem 14.21 states that the Fourier transform 
values are concentrated only on the unit circle f = 1 and converges to zero elsewhere 
as t — y 0, it should be sufficient if the integration for f is done over a region very 
close to r = 1. The following theorem — the principal result in this paper — confirms 
our observation. 

Theorem 5.1. For any given < e < A < S < I, uj G [0, 2ir) and < A < 2ir, 
(5.3) lim / / P^{r,Lo)rdrdu} = / P(u)dw. 



1-6 



5.1. An Intuitive Examination of Theorem 15.11 Before we proceed with 
the formal proof, we again try and give an intuitive explanation of why the theorem 
statement is true. The Fourier transform of the CWR defined in Equation 14.41 involves 
two spatial integrals (over x and y) which are converted into polar coordinate integrals. 
The squared magnitude of the Fourier transform (power spectrum), P^(f, w), involves 
multiplying the Fourier transform with its complex conjugate. The complex conjugate 
is yet another 2D integral which we will perform in polar coordinates. As the gradient 
density function is one-dimensional and defined over the space of orientations, we 
integrate the power spectrum along the radial direction close to the unit circle f = 1 
(as S — >• 0). This is a fifth integral. When we poll the power spectrum P*(r, w) close to 
f = 1, the two sinusoids, namely, exp ^l£i£iEi^ ailc l ex p ^ -i(uoc+vy) ^ ^ Equation 14.41 

are in resonance only when there is a perfect match between the orientation of each 
ray of the distance transform S(x, y) and the angle of the 2D spatial frequency (a; = 
arctan (-))• All the grid locations (x, y) having the same gradient orientation 

(5.4) arctan ( — ^ | = arctan ( — ~\ 

\S X J \u/ 

cast a vote only at their corresponding spatial frequency "histogram" bin us. Since 
the histogram bin is generally populated by votes from multiple grid locations, this 
leads to cross phase factors. Integrating the power spectrum over a small range on the 
orientation (constituting the sixth integral) helps in canceling out these phase factors 
giving us the desired result when we take the limit as r — > 0. This integral and limit 
cannot be exchanged because the phase factors will not otherwise cancel. The proof 
mainly deals with managing these six integrals. 

5.2. Proof of Theorem [5711 We now provide the formal proof of Theorem [5Tj 
For the sake of readability, we divide the proof into smaller subsections. To achieve 
a good flow, we state major portions of our proof as lemmas whose proofs are given 
in the appendix. We would like to emphasize that these lemmas are meaningful only 
within the context of the proof and do not have much significance as stand-alone 
statements. Important symbols used in the proof are adumbrated in Table [5TT1 
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Table 5.1 
Table of important symbols 



Symbol 


Comments 


I 


Integral of P* over the radial length [1 — 6, 1 + 8]. 


9jk 


Integral over the variables r, 9 and f after symmetry breaking. 




Phase term in the integral for gjk- 




Integrals for the main and the error terms of gjk respectively. 
I is the sum of and /■? . 

J— J— . . 


p,q 


Functions used in the definition of . 
p represents the phase term. 


r(ij r(2) 7(3) 
J jk ' J jk ' J jk 


Integrals obtained when is split over the integral range for 9' . 


P 


Symbol used to divide the integral range for 9' into three integrals. 
The limit as /? — > is considered in the proof. 


Gjk 


Result of integrating over 9' while evaluating jj^ . 


tpjk,X 


Integrals for the main and the error terms of Gjk respectively, 
jjjjj is the sum of ipjk and \- 


£3, £4 


Error terms used to define \- 



First, observe that 



(5.5) F T e (f,w) = > — — / exp [l-rcos(6»'-ij)] )r'dr'd9'. 

k T 1 2lTTle Jo Jb£\9') V T 



Define 

pl+S pl+6 

(5.6) I(^)= P*(r,u)rdr= F*{?,uj)F*{?,uj)?d?. 

J is J is 



As r — > 0, we show that I(oj) approaches the density function of the gradients of 
S(x, y). Note that the integral in Equation [STB] is over the interval [1 — 5, 1 + 5], where 
S > can be made arbitrarily small (as r — > 0) and this is due to Theorem 14. 2 1 

Recall that in order to evaluate I(ui), we need to perform five integrals, four 
to obtain the power spectrum P£(r,u)) and a fifth along the radial direction f over 
[1 — 6, 1 + S] which is close to the unit circle f = 1. An easy way to compute I(ui) 
in the limit r — ► would be to apply a 5D stationary phase approximation |23| . 
Unfortunately, the 5D stationary phase approximation cannot be directly employed 
in our case for reasons detailed in Appendix |B| 

Breaking the symmetry of the integral. As described in Section IB. 11 we 
propose to solve for in Equation 15.61 by breaking the symmetry of the integral. 
We fix the conjugate variables r' and 9' and perform the integration only with respect 
to the other three variables namely r, 9 and r. To this end, let 

JL JL 1 r 27r r R i 2) ^ /-j>'\ 
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where 



1+6 p2ir rRf ) (6) 



{5.fy jk (r',6';u>) = I I I exp{~~/ jk (r,8,r;r',6',u>)\ f 2 (r,r)drd6df. 

Ji-s Jo JR^ie) l T J 

Here, 

7jfe(r, 8, f; r', 8' ,v) = r [1 — rcos(8 — w)] + r'rcos(9' — u) 

(5.9) -f [cos(u)(xj - x k ) + sin(w)(y J - y k )] 

and 

(5.10) / 2 (r,r) = rr. 

In the definition of jj k (r,6,f;r' ,8' ,u>) in Equation 15.91 uj,r' and 0' are held fixed. 
Similarly, in the definition of gjk{r' , 9'; ui) in Eguation l5.81 w is a constant. The phase 
term of the quantity CjCk (Equation [B3| is absorbed in jjk and pursuant to Fubini's 
theorem [7] , the integration with respect to f can be considered before the integration 
over r' and 8' . Define 

(5.11) r jk (r',8';uj) = r'cos(6>' - ui) - [cos(u))(xj - x k ) +sin(w)(y J - y k )]. 

This leads to the following lemma. 

Lemma 5.2. If r jk (r' ,8';lu) > 0, then as r — > 0, 

g jk (r',8';u) = (2nr)^r, k (r',9';co)e^ ^ r '^ + ™\ 

(5.12) +r K t; jk (r',8';Lo) 

where n>2 and S,j k (r' ,9' ';w) is some bounded continuous function which includes the 
contributions from the boundary. Ifrj k (r',8';uj) < 0, then as t — > 0, gj k (r', 9'; to) = 0. 



The proof of Lemma 15.21 — obtained using a three dimensional stationary phase 
approximation — is available in Appendix [C] Note that for j — k and 8' close to ui, 
rj k (r\ 9'\ u>) > and hence gj k (r' ,8'; u>) =^ 0. Below, we show that the only pertinent 
scenarios that need consideration are 8' close to uj and j = k. When 8' is away from 
w or j ^ k, the integral gj k (r' , 8'; ui) vanishes. Hence, for the sake of readability of 
our proof, we let gj k {r' ,8' ';w) take the most general form given in Equation 15.121 for 
all values of r' and 8' . 

Determining I(ui). Substituting the value of gj k (r' ,8';uj) into Equation 15. 7\ as 

t — > 0, we get 

(5.i3) ^) = EE{^ 1 ^ ) h+^ ) h| 



where 



j=i k=i 



J \/2ttt Jo Jr^(6') \ t J 

<o, r 27r r< 2} ^ f-ir'\ i 

! > )= L L Pm exp (— y^-^.v-.u)^, 
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f~ia jk (uj) in" 
r) jk (u) = exp I ^ + — 

ajfc(w) = cos(uj)(xj - .T fe ) + sin(a;)(?/j - yu), 
p(r', 9'; w) = -r'[l - cos(0' - w)] and 

(5.14) g(r', 0'; lj) = r'^r 1 cos(6>' - w) - c^w). 

In the definition of the functions p(r',9';ui) and q(r' , 9' ; ui) , uj is held fixed. Since 

(2) 

k > 2, by the Riemann-Lebesgue lemma, we have lim r _>o I)f. = an d from the 
Lebesgue dominated convergence theorem, it follows that 

rWr,.,\ - / h™ tWi 



(5.15) lim / = / lim/^H = 



Using the above result in Equation 15.131 we get 

(5.16) lim I( W )^ = ££]im/ ^/«H^. 

Splitting the integral over 6*' into three disconnected regions. Consider 
the integral jjjjp(cj). As essential contributions to it come only from the stationary 
points of p(r r ,6';v) [HJ OH [3] (with w held fixed), we first determine its critical 
(stationary) point(s). The partial derivatives of p(r',9';u) at a fixed value of w are 
given by 

(5.17) = + |l = _ r 'sin(tf'- W ). 

For Vp = 0, we must have 6' = lo. Hence, in order to evaluate /j^ (u>), we find it useful 
to divide the integral range [0, 2ir) for 9' into three disjoint regions namely [0, lu — (5). 
[ui — (3, u + f3] and (u + /?, 2tt) for a fixed (3 > 0, and write 

(5-18) J#( W ) = J«(/3, W ) + Ji 2) (/3, W ) + J#(/3, w) 



where 



J${P,u) = -±= / exp *™ q(/,9'; U )dT'd0', 

V2TTT Ju-13 Jr£>(0') \ T / 

jgUM . ' f /"''""'exp (Mf^M) mcl 
J V2ttt Jo Jr^($>) \ T / 

(5.19^^,0,) = -^= T r'^^xpfM!^),^,,^)^. 
V2ttt J u +pJrW(0') \ T J 

Since the above relation is true for any /3 > 0, we can let /3 — ?> (after we take the 
limit r — > 0). Fix a (3 close to zero and consider the above integrals as r — > 0. Then 
we obtain: 

Lemma 5.3. 

(5.20) lim f 0+A 2^ig)(a,)cfc, = lim lim J&MjWfo „)*,,. 

The proof is available in Appendix [D] 
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Interchanging the order of integration between r' and 9' . We now evaluate 
jj^(/3,w) by interchanging the order of integration between r' and 9' which requires 
us to rewrite 9' as a function of r'. Recall that the boundaries of each T>\ along 
r{9') = b!u'{Q') and r{9') = R^'(0') respectively are composed of a finite sequence 
of straight line segments. In order to evaluate J^(/3,u), we need to consider these 
boundaries only within the precincts of the angles [ui — /3, oj + 0\ at each Z>|. But 
for sufficiently small j3, we observe that for every ui € [0, 2ir), when we consider these 
boundaries (along R^\o') and R^\9') respectively) within the angles [ui — j3,u> + /3], 
they are composed of at most two line segments as portrayed in Figure 15.11 




Figure 5.1. Boundary considered within the angles [uj — /3, u + /3] is comprised of at most two 
line segments L\ and L2 ■ 



Over each line segment, r'(8') is either strictly monotonic (strictly increases or 
strictly decreases) or has exactly one critical point (strictly decreases, attains a min- 
imum and then strictly increases) as described in Figure \57I\ 




Figure 5.2. Plot of radial length (r) vs angle (9). 



Hence, it follows that for sufficiently small /3, 0' rewritten as a function of r' is 
composed of at most three disconnected regions (as seen in Figure [573| . 

Let B(r') C [ui — f3, lu + f3] denote the integration region for 9'(r'). Treating 9' as 
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Three disconnected regions 
Figure 5.3. Three disconnected regions for the angle (0). 



a function of r', the integral Jj^(0,lo) can be rewritten as 



(5.21) 
where 



Gj k (r',ui)dr', 



r i k 1) (f3,u J )=mi{R { k 1 \e')}, and 
(5.22) rWG8, W ) = B up{i#> (9')} 

with 0' € [w - f),u) + 0\ and 

1 f fip{r',9',uj 



(5.23) 



/27TT 



exp 



B(r') 



q(r',0',u)d9'. 



Note that while evaluating the integral Gjk{r',uj), r' and w are held fixed. As con- 
tributions to Gjk(r',u)) come only from the stationary points of p(r', 9' ,oS) (with r' 
and to held fixed) as r — !• 0, we evaluate 



require 



Moreover 



lo) and for it to vanish, we 



d 2 p 



d9' 2 

p(r' , uj, lu) = 0, and 



(5.24) 



q(r',uj,uj) = r'yjr' - a jk {u)). 



For the given r', if w ^ S(r'), no stationary points exist. Using integration by parts, 
Gjk(r',uj) can be shown to be 63(7-', w,r) = 0(y/r), which can be uniformly bounded 
by a function of r' and oj for small values of r. 

If w € B(r'), then using the one dimensional stationary phase approximation 
[T51 [TB] , it can be shown that 



Gjk(r',w) =exp ( -t^-J Vr'\Jr l - a jk (u}) + e 4 (r', w, t), 



(5.25) 



where 64 (r', w, r) can be uniformly bounded by a function of r' and w for small values 
of t and converges to zero as r — > 0. Here, we assume that the stationary point 9' = lu 
lies in the interior of B{r') and not on the boundary as there can be at most finite 
(actually 2) values of r' (with Lebesgue measure zero) for which 9' = oj can lie on the 
boundary of B(r'). 
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Computing the integral over lo and r' . Let ri \f3,u) > r£\fi,u>) and 
(J3, w) < r[ 2) (j3, lo) be the values of r' such that when (j3, lo) < r' < rjf 1 ^ (/3, lo), 
the stationary point 9' = ui lies in the interior of B(r'). Substituting the value of 
Gjk(r',ui) into Equation 15.211 and using the definitions of r]jk(w) and ajk(Lo) from 
Equation 15. 14i we get 

(5.26) 
where 

1 ^^o+a /-i a .,,(u)\ r r<k+){J3 ^ ) ,— i 

(5.27) Vjfc(£) = 77 / exp / vWr' - a ]k (uj)dr'dw 



and 



e 4 (r>,T), r<f } (/3,w) < r' < r{ +) (/3,w), 
rj^GS.w) or r{ 4 



1 eafr'.w.T), r' < rj^GS, lo) or 4 +) (/3, w) < r'. 



Since ^^(a;)! = 1 and x(r',w, r) can be uniformly bounded by a function r' and 
lo for small values of r, by the Lebesgue dominated convergence theorem we have 

lim n +A I!ilMl [ rk iM ' x (j',u,T)dr\duj 
(5.28) = / { / lim x{r',cJ, r) dr' \ duo = 0. 



This leaves us having to prove the following result: 
Lemma 5.4. 

(5.29) Y, E ] im ]To ^ifcCS) = / i'Mdw- 

3=1 fc=l T ^ •'"D 

The proof of this lemma is given in Appendix [E] This completes the proof of Theo- 
rem 15.11 

We would like to give a short recap of our proof. Beginning with the defini- 
tion of I(lo) in Equation 15.61 Lemma \5. 21 and the statements following it lead to the 
relation 15. 161 namely, 



-in j =1 k=1 T Juin 



From Lemma [5751 it follows that 

K K r u +A 



E E ]To / ~~L^ fc {uj)du = 

K K /•wo+A / \ 

(5.31) EE lim lim / ^2 jj^, w)du . 



j=i fc=i 



1G 



KARTHIK S. GURUMOORTHY AND ANAND RANGARAJAN 



Interchanging the order of integration between r and 9' , we showed that 

K K ruo+A r \ K K 

(5.32) EE" ^ J f^M^ = Y,Y,}™ H^- 

j=lk=i ua i=ifc=i 

Finally, the application of Lemma |5 . 41 gives the desired result of Theorem 15. 1[ namely, 

r i*j +A pl+S /■wo+A 



lim / / P^(r,uj)rdrd(jj = lim / I(u>)duj 

^ cun J x — o J ujq 

K K 

= V" V" lim lim ipj k {/3) 



i=X fe=l' 

(5.33) = / P{u>)du>. 



6. Results Stemming from the Main Theorem. As an implication of The- 
orem [3TH we have the following corollary. 

COROLLARY 6.1. For any given < 6 < 1, w € [0, 27r), 



I /-cjo+A f f.l+8 1 

(6.1) lim lim — lim / < / P*{r,uS)r df > duj = P(ojo). 

£ ->0A-+(1At->I)J Uo I I 



Proof. From Equation 13.51 we have 

(6.2) lim I P +A P (w) cL, = hm ^^i^) = P(wo) . 
A^oAy wo a^o a 

Since Theorem 15.11 is true for any < e < i, it also holds as e — > 0. The result then 

follows immediately. 

□ 

Theorem 15.11 also entails the following lemma. 
Lemma 6.2. For any given 0<e<|, 0<<S<1, 

(6.3) lim/ / P*(r,u))rdrdw = 1. 

Proof. Since the result shown in Theorem 15.11 holds good for any loq and A, we 
may choose u>q = and A = 2ir. Using Equation 13.61 the result follows immediately 

as 

(6.4) lim/ / P^(f,uj)?dfduj^ P(bj)duj = 1. 
□ 

Lemmas 16.21 and 14.11 leads to the following corollaries. 
COROLLARY 6.3. For any given 



(6.5) 



/>27r f /"l— 8 poo j 

lim y j y P T £ (f , w) fdf + J P T £ (f , w) fdf > do; = 0. 
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Proof. From Lemma l4.11 we have for any r > and < e < |, 

(6.6) / / P*(f,u))rdrdu) = 1. 

./o 

For the given < 5 < 1, dividing the integral range (0, oo) for r into three disjoint 
regions namely (0, 1 — 6), [1 — 6, 1 + 6] and (1 + 6, oo) and letting r — > 0, we have 

/ 1 27T | />1 — 5 />00 | 

lim / < / P^(r,u) fdf + / P^(r,cv) fdf + / P f f (r»fdf Ui=l. 



Pursuant to Lemma |6.2[ the limit 



(6.7) lim / / Pr(f,u) fdfduj 



T-S-O 



2tt /.1+5 



J 1-8 



exists and equals 1. The result then follows. 
□ 

COROLLARY 6.4. For any given < e < |, < 6 < 1, Wo € [0, 27r) and 
< A < 2tt, 



T->0 







lui 


L 



1-5 



(6.8) lim/ < / P*(r,u)rdr+ P^w) fdf } du; = 0. 



1+5 



Proo/. Let M = [^J . Define uj i+1 = lo, + A mod 2tt for < i < M - 1. Then, 
we have from Corollary [ 



(6.9) lim 

T— >0 



where 



M— 1 «u;i-)-i />a;o+27r 

^ / Q(w)dw + / Q{uj)duj 



= 0, 



/■l — <5 /-DC 

(6.10) Q(uj)= P e T {f,uj)fdf+ l P T e (f , w) fdf . 

io il+5 

Since P^(f,Lu)f > 0, it follows that Q(w) and both integrals in Equation l6.9l are non- 
negative and hence each integral converges to zero independently giving us the desired 
result. 
□ 

From Theorem 15.11 and Corollaries 16.11 and 16. 4\ the subsequent results follow almost 
immediately. 

PROPOSITION 6.5. For any given < e < \, cvo S [0, 2ir) and < A < 2ir, 

(6.11) lim / \ Pr(r,u)r dr\ dw = / P(uj)duj. 



COROLLARY 6.6. For any given w G [0, 2ir), 

p^o+A ( f oo 



.12) lim lim — lim / 11 P?(f,u)r df\duj = P 



(w ). 
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7. Significance of our result and concluding remarks. The integrals 

(7.1) / / P^(f,uj)Mrdu}, / P{uj)duj 

give the interval measures of the density functions P^ (when polled close to the unit 
circle f — 1) and P respectively. Theorem 15 . 1 1 states that at small values of r, both the 
interval measures are approximately equal, with the difference between them being 
o(l). Furthermore the result is also true as e — > 0. Recall that by definition, P^ is 

the normalized power spectrum of the wave function 4>{x, y) — exp ( iS &< v ^ . Hence, 

we conclude that the power spectrum of (j)(x, y) when polled close to the unit circle 
f = 1 (as 5 — > in Theorem 15.11) , or when integrated over r (with reference to 
Proposition 16 . 5p . can potentially serve as a density estimator of the orientation of V5* 
for small values of r and e. Our work is essentially an application of the higher-order 
stationary phase approximation culminating in a new density estimator. 

7.1. Advantages of our formulation. One of the foremost advantages of our 
method is that the orientation gradient density is computed without actually deter- 
mining the distance transform gradients. Since the stationary points (as seen from 
the stationary phase approximation) capture gradient information and slot them into 
the corresponding frequency bins, we can directly work with the distance function — 
circumventing the need to compute its derivatives. We are not aware of any previous 
work that estimates the orientation gradient density without first computing the gra- 
dients of the distance transform. 

Recall that we furnished a closed-form expression for the distance transform gra- 
dient density function P{9) in Equation 13.51 While it initially appears attractive, 
computing the density function via the closed-form expression is practically cumber- 
some as we need to first determine the Voronoi region corresponding to each Voronoi 
center (source point) Yfc and then for each orientation direction 0. compute the ray 
length Rk{9) from Yfc to its Voronoi boundary along 0. These involve unwieldy ma- 
nipulations of complex data structures. On the other hand, our mathematical result 
provides an easy mechanism to achieve the same task as it is computationally faster 
and easier to implement. Given the N sampled values S of the distance function S 
from a point-set of cardinality K, we just need to compute the fast Fourier trans- 
form of exp ( ^fe) j — an 0(N log N) operation — and then subsequently compute the 

squared magnitude (to obtain the power spectrum) — performed in O(N). Hence the 
orientation density function can be determined in 0{N log N) independent of the car- 
dinality of the point-set (K). Our algorithm is computationally efficient even when 
K = O(N). 

7.2. Possible Extensions. The present work only deals with special kinds of 
distance functions, namely those defined from a set of discrete point locations in two 
dimensions. Other cases include signed distance functions, distance functions defined 
from a set of curves |18l |2"0] etc. While our work initially appears to be somewhat 
restrictive, we should note that as the cardinality and locations in each point-set can be 
arbitrary, the resulting distance functions can be quite complex. We strongly believe 
that it is possible to establish a similar Fourier transform-based density estimation 
result even for arbitrary, continuous (and differentiable) functions in 2D. A general 
result of this nature would subsume our current work on point-set distance functions 
as well as the other kinds of distance functions mentioned above, as the gradient 
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magnitude of an arbitrary 2D function can vary across the point locations and need 
not necessarily be identically equal to one as in the case of distance functions. Here 
the gradient density function will inherently be two dimensional and defined both 
along the radial, gradient magnitude direction (r = \J + S^j and the orientation 

^9 = arctan )J • However, at the present time, due to many technical issues, this 
is merely a conjecture and requires further investigation. Generalization of the present 
work to three dimensions is also a concrete possibility. These are fruitful avenues for 
future research and we may explore them in the years to come. 

Appendix A. Proof of Lemma 14.11 

Proof. Define a function H(x,y) by 



Let 



H(xv) = ( 1: ^ (x ' 2/)e " £; 
y ,yj [ : otherwise. 

f(x, y) = H(x, y) exp (i^al) . Then, 



(A.l) F<( U) v) = -^- e JJ /(*,!/) e*p ( + dxdy . 

Let 2 = Sj V. -t an d G(s,t) = F*(sT,tr). Then, 

(A.2) tZ £ G(s, t) = -^ff f(x, y) exp (-i(sx + ty)) dxdy. 

Since / is t 1 integrable, by Parseval's theorem [1], we have 

(A.3) JJ \f(x,y)\ 2 dxdy = jj \rl e G(s,t)\ 2 dsdt = (W e ) 2 jj |F r e (sr, ir)| 2 dsdt. 

Letting u = sr, v = tr and observing that 



(A.4) jj \f(x,y)\ 2 dxdy = // |exp (^MTj 

we get 

(A.5) (Z e ) 2 II '\F*(u,v)f dudv = L e 



2 



dxdy = L c , 



Hence 



2 dudv 



(A.6) jj \F*{u,v)\ 

which completes the proof. 
□ 

Appendix B. Difficulty with the 5D stationary phase approximation. 

Since P*(r,cd) equals F*(r,u))F*(r,u>), we have 

( B -i) f ( w )=EE»w 

j=i k=i y ' 
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where 

(B.2) N jk (u) = 
Here. 



1+6 r 2ir r R { k '(e') r 2n r R { k '(0) 



1-5 JO 



R^He') Jo 



exp ( — bjk I fidrdOdr' ' dQ 1 dr. 



b jk (r,O,r',0',r;u)) = r [1 - f cos(0 - w)] - r' [1 - f cos(0' - lj)] 
(B.3) -f [cos(cj) {xj-x k ) + sin(cj) - y fe )] 

and 

(B.4) h{r,r',f) = rr'f. 

Notice that the phase term of the quantity Cj C'k , namely 
(B.5) - r [cos(a;)(a;j - x k ) + sia(w)(% - y k )] 

is absorbed in bj k . Since we are interested only in the limit as r — > 0, essential 
contribution to Njk(ui) comes only from the stationary (critical) point(s) of bjk |23) . 
The partial derivatives of bj k (r, 6, r',9', r) are given by 



db 



(B.6) 



dr 
dbjh 
dr' 
db jk 

dr 



1 — fcos(9 — oj), 
-1 + f cos(0' - lj) 



db 



3k 



89 
db jk 
89' 



rf sin(0 — lj), 

= —r'fsialO' — u>), and 



-r cos 



w) + ?*'cos(0' — lu) — [cos(ui)(xj — Xk) + sin(w)(?/j — Vk)]- 



As r, r' and f > 0, it is easy to see that for Vbj k = (stationary), we must have 

(B.7) f=l, = 6' = ui,r = r' - [cos(uj)(xj - x k ) + sin(cj)(^ - y k )]. 

Let to denote the stationary point. The Hessian matrix VV of bj k at to is given by 



W(r,9,r',e',f)\ t 



where r to = r' — [cos(oj)(xj — x k ) + sin(u;)(t/j — y/c)]. Unfortunately, the determinant 
of W at the stationary point io equals as the first and third rows — corresponding 
to r and r' respectively — are scalar multiples of each other. This impedes us from 
directly applying the 5D stationary phase approximation [23J. 

The addition of a 6 th integral to the above setup — where the power spectrum is 
integrated over a small range on the orientation u> (in order to remove cross phase 
factors) — leads to a 6D stationary phase approximation. This is of no help either, as 
the Hessian continues to remain degenerate for the same reasons as above. 
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B.l. Avoiding degeneracy by symmetry breaking. As we notice above, 
the degeneracy in the 5D (and 6D) stationary phase approximation arises because 
the determinant of the Hessian, namely W, when evaluated at the stationary point to 
takes the value zero, as its first and third rows corresponding to r and r' respectively 
are scalar multiples of each other. Also, observe that the value of either r or r' is not 
determined at the stationary point and can take on arbitrary values. However, the 
rows (and columns) of W corresponding to the other three variables 9, 9' and f are 
indeed independent of each other and do not cause degeneracy. This strongly suggests 
that if we do not consider both r and r' together and hold back either one of them, 
say r' , the resulting 4D stationary phase approximation will be well-defined. Since 
the integration range for r' is defined in terms of 9' , we retain both these variables 
and perform the stationary phase approximation on the other three variables. This 
manual breaking of symmetry avoids the degeneracy issue. 

Appendix C. Proof of Lemma 15.21 

Proof. Recall that the essential contribution to gj k (r',9';cj) comes only from the 
stationary points of "fj k as r — > |23| . The partial derivatives of jj k (r,9,r;r J ,9',u)) 
are given by 

1 — rcos(9 — lj), J — rrsui(9 — oj), and 
o9 

—rcos(9 — u) + r' cos(0' — lj) — [cos(ui)(xj — x k ) + sin(cj)(yj — y k )]. 

As both r and f > 0, for V^jk — (stationary), we must have 

f = 1, 9 = ijJ, and 
(C.2) r = r' cos(0' - lj) - [cos(w)(a;j - x k ) + sin(o;)(yj - y k )]. 

Let to denote a stationary point. Then 

7ifc(*o) = r' cos(9' - lj) - [cos(w)(x J - x k ) + sm(u))(yj - y k )\ = r jk (r' , 9'; lj), 
/a(*o) =r jk (r',9';oj) 

and the Hessian matrix % of at the stationary point to is 

-1 " 

r jk (r',9';Lj) . 
-10 

It can be easily verified that the determinant of H equals —rj k (r', 9'). 

If Tj k (r', 9'; lj) < 0, no stationary points exist as r > by definition and hence 
gj k (r' ,9';lj) — as r — > [23J. If rj k (r' ,9';lj) > 0, the determinant of % is strictly 
negative and its signature-the difference between the number of positive and negative 
eigenvalues-is 1. Then, from the higher-order stationary phase approximation [23J, 
we have 

gMr',0';Lj) = (2nr)i yV^r', 9'; lj) exp + + ^ j ^ 

as r — > 0, where t\(r' , 9' , r; lj) includes the contributions from the boundary in Equa- 
tion 15.81 Here we have assumed that the stationary point to does not occur on the 



(C.l 



dr 
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boundary and lies to its interior, i.e, R^ (9) < rjk{r' ,6'\w) < Rj(6), as the measure 
on the set of {uj,9',r'} for which rjk(r' ,9';uj) (or to) can occur on the boundary is 
zero. 

Let r denote the boundary in Equation l5.8l If there does not exist a 2D patch on 
r on which "fjk is constant, then we can conclude that ei(r', 9', t; uj) — which includes 
the contributions from the boundary T involving the stationary points of the second 
kind where the level curves of are tangential to T — should be at least 0(t 2 ) as 
t -)• [Hid [53]. From this, we get 

(C.3) ei(r',6',T;u,)=T% k (r',6';oj) 

where k > 2 and ^jk(r' , 9'; uj) is some bounded, continuous function. Since the bound- 
ary r is made of straight line segments, we can show that this is indeed the case. 
Below, we take care of this technical issue. 

The boundary T in Equation 15.81 is the union of two disconnected surfaces Ti = 
Ai x [1 — S, 1 + 6] and T 2 = A 2 x [1 — S, 1 + S] where A\ is the boundary along 
r(9) = Rf\0) and A 2 is the boundary along r(9) = Rf\6). Note that both Ai and 
A 2 are composed of a finite sequence of straight line segments. Consider the surface 
Y\. The value of jj k on the surface Ti at a given 9 and f (with r' , 9' and uj held 
fixed) equals 

(C.4) t£ (9, f; r', 9', uj) = pf ] (9) [1 - f cos(9 - uj)} + rr jk (r\ 9';uj). 

Following the lines of Theorem 14. 2 [ we observe that for a given f, 7^ (9, f; r' , 9', lu) 
cannot be constant for a contiguous interval of 9 as Equation 14. 181 cannot be satisfied 
over any finite interval. By a similar argument, there can exist at most only a finite 
discrete set of 9 for which R i 'j l \9) cos(# — uj) — rjk(r' , 9';lu). Let Z denote this finite 

set. Then, for a given 9 ^ Z, jj£ varies linearly in f and specifically, its derivative 
with respect to f does not vanish. From the above observations, we can conclude 
that there does not exist a 2D patch on Ti on which 7^ is constant. A similar 
conclusion can be obtained even for the surface Hence, 7^ cannot be constant 
on the boundary T over a 2D region having a finite non-zero measure. 
□ 



Appendix D. Proof of Lemma [ 

f2l (3) 

Proof. By construction, the integrals J^ k (/3, a;) and J. fc (/3, lu) do not include 
the stationary point 9' = uj and hence Vp 7^ in these integrals. Following the 
lines of Theorem 14. 2 [ by defining the vector field u = ir^jp g and then applying the 

divergence theorem, both jj^(/3,w) and J^\/3,uj) can be shown to be r E2 (' 2 ' (/?, lu) 

and t K3 (^ 3 ) (/?, w) respectively where both k 2 and K3 > 0.5 and <^ 2 ) and C*- 3 -* are some 
continuous bounded functions of /3 and lu. Hence, we can conclude that 

r 2ir 



(D.l) 



>0 







,. T 

< lim — 

t->o L 



><■> 







C (2) (A-) 



duu = 



( 31 

as \r)jk = 1| and similarly for J. fc (/3, u>) for any fixed /3 > 0. It follows that the result 
also holds as j3 —> provided the limit for /3 is considered after the limit for r, i.e, 

^J J ( 2 ) G8,w)du = Q, and 
(D.2) lul !, ,nl , l 1 / ^-^(M^^. 







lim lim 












lim lim 






lu 
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Hence, I^(uS) in Equation 15.181 can be approximated by J^(P,lo) as j3 — > and as 

r ->• 0. 

□ 

Appendix E. Proof of Lemma 15.41 

Proof. Define 



(E.l) p jk (p,w)= Vr'Jr' -a jk {u)dr'. 

We consider two cases, one in which j — k and another in which j =/= k. 

case(i): If j ^ k, then aj k (ui) varies continuously with lo. Also, notice that 
Pjk{P, w) is independent of r and is also a bounded function of /3 and lo. The stationary 
point(s) of ajk — denoted by Co — satisfy 

(E.2) tan(w) = Vi ~ Vk , 

and the second derivative of ayfe(cj) at its stationary point(s) is given by 
(E.3) a%(Co) = -a jk (Q). 

For a'j k (u)) — 0, we must have 

(E4) tan(w) = -El^h. = ^JL^ , 

% - yfc atj - at* 

where the last equality is obtained using Equation IE.2I Rewriting, we get 

(E5) c^i* y = _i 

which cannot be true. Since the second derivative cannot vanish at the stationary 
point Co, from the one-dimensional stationary phase approximation [T5] , we have 

(E.6) lim — r° +A exp ( ~ ia 3 k ^ \ Pjk (B u )du = lim 0(t k ) = 

L e J uo V T J 

where k — 0.5 or 1 depending upon whether the interval [loq,loq + A) contains the 
stationary point (Co) or not. Hence, we have ipj k ((3) — for j k. 
case(ii): If j — k, then a kk [ui) = and 

r< +) (/3,^) 



p kk (p,u>) = r dr , 



(E.7) ipkk{0) = jj I pkk(/3,Lo)duj. 



U-'O 



From the definitions of r^\p,Lo) and r^(P,tS) in Equation 15.221 we observe that 

limr^&wJt^M, and 
(E.8) lim rf^H^V)- 
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Since r { ~\l3,Lu) -> r^\/3,ui) and r ( k +) {/3,uj) -> rf\p,w) as /3 -> 0, we have 

limr^^H^V), and 
(E.9) limr< +) (/J,a;)=fl<%). 

Since r^ } (/3,w) > r[ 1} (/3,w) and r k + \f3,uj) < r^ 2) (/3,o;) at a fixed (3 and r' > 0, we 
see that Pkk(P, w) can be bounded from above by a positive decreasing function of /?, 
namely, 



(E.10) PkkGS.w) < /" ' 



r'dr', 



and is also independent of r. As both r^(/3, w) and (/3, uj) are also bounded 
functions, by the Lebesgue dominated convergence theorem, we get 

lim lim Vfefe(/3) = — / Hm p fefe (/3,w)dw 

"^Lo {Jul 1 

(E.n, _ &i *>i»* ! S£ 



r'. 



Recall that L e = (1 - 2e)£. Hence 



j=ik=i k=i Ju) ° 
(E.12) = / P(w)dcj 



which completes the proof. 
□ 
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